Oscillatory zoning of minerals as a fingerprint of impurity-mediated growth

We propose a kinetic mathematical model of the oscillatory compositional zoning profile recorded in minerals based on the crystal growth suppression induced by impurities. Notably, the presence of a small amount of impurities significantly inhibits crystal growth, and a growth inhibition mechanism called the pinning effect is widely accepted. Here we show that a model that considers the pinning effect and adsorption/desorption kinetics of impurities on the crystal surface can reproduce the oscillatory compositional zoning. As impurities are common in nature, this model suggests the existence of a universal mechanism that can occur in the growth processes of various crystals.


Oscillatory zoning of minerals as a fingerprint of impurity-mediated growth Hiroki Torii & Hitoshi Miura *
We propose a kinetic mathematical model of the oscillatory compositional zoning profile recorded in minerals based on the crystal growth suppression induced by impurities.Notably, the presence of a small amount of impurities significantly inhibits crystal growth, and a growth inhibition mechanism called the pinning effect is widely accepted.Here we show that a model that considers the pinning effect and adsorption/desorption kinetics of impurities on the crystal surface can reproduce the oscillatory compositional zoning.As impurities are common in nature, this model suggests the existence of a universal mechanism that can occur in the growth processes of various crystals.
A zonal structure is observed on mineral sections wherein different compositions show nearly parallel distributions.Among these structures, those wherein compositional changes are periodically repeated are called oscillatory compositional zoning (OCZ).Plagioclase, a major rock-forming mineral, is a notable example exhibiting OCZ 1 .Plagioclase is a solid solution with anorthite (An; CaAl 2 Si 2 O 8 ) and albite (Ab; NaAlSi 3 O 8 ) as end members, and the mole fraction of An varies periodically.OCZs have also been identified in Ti in quartz 2 and trace elements, such as P and Al, in olivine 3 .OCZs have also been observed in pyroxene for major and trace elements 4 .Elucidating the conditions and mechanisms of OCZ formation is important for deciphering the formation history of these minerals.
The cause of OCZ can be classified into two main categories: external (periodic fluctuations of external factors, such as temperature, composition, and pressure) and intrinsic (mechanisms inherent in the crystal growth process itself).Here, we focus on the latter, especially on the crystal growth in solution.In the intrinsic mechanism, the influence of a compositional boundary layer formed around the growing mineral is considered 5,6 .The elemental partitioning between the mineral and ambient solution causes compositional changes in the boundary layer and subsequently changes the driving force of crystallization at the mineral surface.The formation of OCZ requires a nonlinear relationship between the driving force and crystal growth rate.Three main models have been proposed to explain this nonlinear relationship.The first model is called surface reconstruction, wherein the crystal growth rate depends on the microstructure of the crystal surface, but the microstructure does not respond immediately to changes in the driving force [7][8][9] .The second model is mainly applied to plagioclase, wherein the rate of incorporation of An or Ab into plagioclase depends on the mole fraction of each component on the surface 10,11 .The third model considers non-equilibrium elemental partitioning at the crystal surface and is based on the assumption that the partition coefficients are reversed from those at equilibrium [12][13][14] .
In this study, we focus on the effect of impurities on crystal growth, which has not been considered in the aforementioned models, and theoretically show that the general crystal growth inhibition effect induced by impurities causes OCZ.Here, impurities are defined as components other than the major crystal constituents present in the solution.

Crystal growth hysteresis induced by impurities
Crystal surface is stacked in layers by atomic steps moving forward, generally generated by screw dislocations exposed to the surface 15 .Impurities are adsorbed on the growing crystal surface and inhibit its growth by checking the advancement of these steps (Fig. 1A).This is called the step pinning effect and has been widely confirmed in the growth process of various crystals 16 .Conversely, the impurity adsorption rate is known to be finite, which causes complex nonlinear phenomena owing to competition with the crystal growth rate.One theory that explains this nonlinear phenomenon is the bistability of the growth rate [17][18][19][20] .Figure 1B shows the dependence of the growth rate on the degree of supersaturation (see Methods for details).At high supersaturation, the crystal surface is kept renewed with new successive steps before it is covered with impurities; thus, few adsorbed impurities are present and the crystal growth is hardly hindered (point a).As the supersaturation decreases, the frequency of surface renewal decreases, adsorbed impurities gradually increase, and the growth rate decreases (path b).When the supersaturation falls below a certain critical value, the step advancement stops completely owing to positive feedback as the surface that is no longer renewed is further contaminated by impurities (point c).This surface is saturated by adsorbed impurities (a state where adsorption and desorption of impurities are balanced, namely, adsorption/desorption equilibrium); thus, the crystal growth does not resume immediately even if supersaturation begins to increase (path d).However, once the growth resumes, the impurity-free surface is restored (point a).The property wherein multiple states can be realized in the intermediate supersaturation region is called bistability, and it explains the hysteresis phenomenon wherein different growth histories are followed by an increase or decrease in supersaturation.

Oscillataory zoning
We conducted simple numerical calculations considering elemental diffusion in the solution to show that the impurity-induced crystal growth hysteresis causes OCZ.Assuming that the crystal surface is an infinitely wide plane, the time-related variation of the concentration distribution of solute and impurity molecules in the solution is obtained by solving the diffusion equation, and the relationship between the interfacial supersaturation and crystal growth rate is elucidated based on the growth hysteresis theory (see Methods for details).The time-related variation of the normal growth rate (the crystal growth rate perpendicular to the crystal surface) obtained from the numerical calculation is shown in Fig. 2A,B.Notably, the growth rate periodically increases and decreases.Figure 2C,D show the time-related variation of the concentration distribution in the solution during crystal growth and growth stoppage, respectively.Notably, a compositional boundary layer is developed near the crystal surface; the concentrations of both the solute and impurity decrease during growth, and are then recovered by diffusive transport from the bulk when the growth is terminated.
Based on above results, we calculated the compositional zoning profile in the grown crystal.The impurity fraction with respect to the solute incorporated into the as-grown surface is proportional to the impurity fraction in the solution.In this study, the impurity fraction in the solution near the crystal surface was multiplied by the constant K to obtain that in the grown crystal at that time (see Methods for details).K is a quantity commonly referred to as the partition coefficient.The resultant zoning profile is shown in Fig. 3 with the blue line.Immediately after the start of crystal growth, the impurity fraction decreases with crystal growth.This is because the impurity diffusion in the solution is assumed to be slower than the solute diffusion in this calculation, so the impurities in the compositional boundary layer are more depleted with respect to the solute molecules.When the crystal growth is then stopped, the impurity fraction increases as the impurities depleted in the boundary layer are supplied by diffusive transport.Thus, the impurity fraction in the crystal at the time of re-growth is higher than that immediately before the termination of growth.Therefore, the impurity fraction in the crystal changes periodically; namely, OCZ is developed.
The characteristics of OCZ, such as interval (wavelength) and amplitude of compositional fluctuations, depend on the impurity adsorption kinetics.For comparison, we demonstrate the calculation for the case where the density of adsorbed impurities at adsorption/desorption equilibrium is the same, but when a longer time is required to reach equilibrium (orange line in Fig. 3).OCZ was reproduced but the wavelength was longer than that of the case when the impurity adsorption is faster.In the case of slow adsorption, the crystal surface is less likely to be contaminated by impurities, so its growth can be maintained until the supersaturation becomes The crystal surface has a step, which is as high as one atom, and the continuous incorporation of atoms into this step moves the step forward and builds the surface layer by layer (layer growth).The step is generated by screw dislocations exposed to the surface.Adsorbed impurities inhibit the step advancement by pinning it in place (pinning effect).(B) Steady-state solution of the step velocity when considering impurities that are repeatedly adsorbed and desorbed on the surface.The abscissa and ordinate represent the degree of supersaturation and dimensionless step velocity, respectively.In an intermediate supersaturation region (bistability region), there exist two solutions: growth (purple) and no-growth (blue) states.The dashed curve indicates an unstable solution that does not appear in practice.lower.Therefore, the wavelength is increased as the crystal can grow for a longer time before being terminated.This indicates that the wavelength of OCZ not only depends upon the diffusion coefficient in the solution but also on the behavior of elements on the crystal surface, suggesting that more careful consideration is required in interpreting the formation condition of OCZs.
The present study shows that OCZs can occur under typical solution growth conditions (see Methods for details).We present qualitative predictions for cases where OCZs do not occur.The first case is when crystal  www.nature.com/scientificreports/growth is not controlled by diffusion.In our model, growth oscillation is triggered by local changes in solute concentration in the boundary layer near the crystal surface.Therefore, the proposed mechanism is ineffective in systems where (1) uptake at the crystal surface determines the growth rate; (2) the free space between adjacent crystals is sufficiently small compared to the diffusion length; or (3) the boundary layer is disturbed by convection.The second case is when the bistability region is sufficiently narrow.Our model assumes the existence of the bistability region (Fig. 1B), and predicts its width (see Methods for details).The bistability region is expected to be very narrow in systems where (1) crystals are made of poorly soluble materials; (2) crystal growth is controlled by the surface incorporation process; (3) the desorption/adsorption of impurities is extremely rapid; (4) the amount of impurities adsorbed on the crystal surface is limited; or (5) the adsorbed impurities are incompletely removed by a single step passage.In such systems, the proposed mechanism will be ineffective as the bistability region only appears under limited supersaturation conditions.The desorption/adsorption behavior of impurities on crystal surfaces is often poorly understood, and it is difficult to compare our model directly with experimental data, however the qualitative predictions discussed herein is usable to validate our model.

Conclusions and outlook
As many impurities are present in nature, the impurity-induced growth hysteresis can occur in general.OCZs found in minerals, such as quartz, olivine, and pyroxene [2][3][4] , may be the result of incompatible elements that are not easily incorporated into minerals playing the role of impurities.The periodic fluctuation of An content in plagioclase 1 may be ascribed to the growth-inhibiting effect of impurities, rather than to the previously proposed growth kinetics unique to binary systems 10,11 .In addition to those of these minerals, the oscillatory behaviours and the periodic structures have been reported in various crystal growth systems that may be attributed to the action of impurities.The growth hysteresis and spontaneous oscillatory growth observed in ice crystal growth experiments are also thought to be related to the effect of proteins added as impurities in the system [21][22][23][24] .Recently, calcium-binding proteins have been found to be periodically distributed in human kidney stones, which are mainly composed of calcium oxalate, and the effect of these proteins on the stone growth has garnered considerable attention [25][26][27] .The impurity effects on crystal growth discussed in this study may be a common cause of the oscillatory phenomena and periodic structures observed in various crystals.

Bistability theory of crystal growth hysteresis
The bistability of the step velocity in the presence of impurities adsorbed on the crystal surface has been studied theoretically 17,18 and numerically 19,20 .Figure 1A shows the schematic of a crystal surface undergoing spiral growth.The spiral step, centered on a screw dislocation, moves forward, building the crystal surface layer by layer.The step interval is given by = 19ρ c , where ρ c = s m κ k B Tσ is the radius of a two-dimensional critical nucleus 28 .s m is the area of the crystal surface occupied by a growth unit, κ is the step-edge free energy, k B is the Boltz- mann constant, T is the absolute temperature, and σ is the supersaturation.If the step velocity is V i , the crystal surface is updated on average with the period τ update = /V i .However, impurities in the solution are adsorbed onto the crystal surface at a certain rate.A faster adsorption process causes the crystal surface to be constantly contaminated (adsorption/desorption equilibrium).Conversely, if the adsorption is slow, the crystal surface is successively renewed before impurities are adsorbed again, thus maintaining a clean surface.This suggests that the density of adsorbed impurities N i (number per unit area) depends on V i and .At the same time, V i depends on N i as the steps are prevented from advancing by the adsorbed impurities.The interdependence between V i and N i clarifies the emergence of crystal growth hysteresis.Herein we outline the theory of growth hysteresis based on the Miura formulation 20 .
The step moves forward by slipping through the gaps between the adsorbed impurity molecules.During this, the steps are bent and velocity is reduced.This mechanism of crystal growth inhibition is widely accepted as the pinning mechanism 16 .The average step velocity V i is as follows 29 .
where V ∞ is the velocity of a straight step when the effect of impurities is negligible and proportional to super- saturation, expressed as V ∞ = β st v m c e σ .β st is the step kinetic coefficient, v m is the molar volume of the solute molecules, and c e is the solubility 30 .We use a coefficient V 0 = β st v m c e to denote V ∞ = V 0 σ .The interior of the square root becomes negative when N i > 1/(2ρ c ) 2 , whereby step advancement is completely inhibited ( V i = 0).
When the step passes, the crystal surface is updated and N i is reset to zero.Subsequently, the adsorbed impu- rities increase with time.The adsorbed impurity density when another step passes is expressed as a solution of the time-dependent Langmuir adsorption equation as follows.
where N e is the value of N i at the adsorption/desorption equilibrium and τ ad is the time needed to reach the adsorption/desorption equilibrium (adsorption time).If the update frequency is high ( τ update /τ ad ≪ 1 ), we obtain N i /N e → 0 , indicating that the crystal surface is updated before reaching the adsorption/desorption equi- librium.If the update frequency is low or crystal growth has stopped ( τ update /τ ad ≫ 1 ), we obtain N i /N e → 1 , indicating that the crystal surface is contaminated with impurities until the adsorption/desorption equilibrium is reached.
Equations ( 1) and ( 2) formulate the interdependence between V i and N i .As the crystal steadily grows, these two equations are satisfied, this is called the steady-state solution.To identify the factors that determine the (1) www.nature.com/scientificreports/steady-state solution, we make these equations dimensionless.The step velocity is normalized as Vi = V i /( 0 /τ ad ) and V0 = V 0 /( 0 /τ ad ) , where 0 is the step interval when σ = 1 and is defined by = 0 /σ .The density of the adsorbed impurities is normalized as θ = N i /N e .Using these normalized quantities, Eqs. (1)and ( 2) are rewrit- ten as follows.
where κ = 2s m N 1/2 e κ k B T is the normalized step-edge free energy.θ is the value obtained by time averaging Eq. ( 2) from immediately after the step pass to the next step pass.Vi = 0 when θ = 1 and σ < κ (Eq.3).Thus, κ corresponds to the minimum supersaturation required for a crystal surface that has stopped growing to start growing again.The relationship between Vi and θ for the three cases with different supersaturations is shown in Supplementary Figure 1.When σ = 0.6 , Eqs. ( 3) and ( 4) have an intersection a with Vi = 0 , indicating that step advancement is completely inhibited (no-growth solution).When σ = 1.1 , there is an intersection b with Vi > 0 , indicating that the step continues to advance (growth solution).Multiple intersections exist, a, c, and d, when σ = 0.9 .Thus, a solution is chosen according to previous growth history.By numerically solving Eqs. ( 3) and ( 4) for the two given dimensionless quantities V0 and κ , we can obtain Vi in the steady growth solution.
Figure 1B shows the calculation results for the case with V0 = 5 and κ = 1 .See the next section for the detailed calculation method.Only the growth solution exists when σ > κ , and only the no-growth solution exists when σ < 0.716 .When 0.716 < σ < κ , both solutions exist (bistability region), thus, the preferred solution depends on previous growth history.

Calculation of steady solution
The point where the growth solution disappears is called the transition point, and the supersaturation and step velocity at this point are σ * and V * , respectively.If p = σ * / κ , q = V * σ * , and u = V * / V0 σ * , then the following holds 20 .Substituting Eq. ( 7) into Eqs.( 5) and (6) to eliminate p, we obtain the following.
The step velocity Vi has multiple solutions when σ * < σ < κ (Fig. 1B).The preferred steady-state solution depends on previous growth history.If the crystal is growing, a growth solution is realized when σ ≥ σ * , and a no-growth solution when σ < σ * .However, when crystal growth is halted, a no-growth solution is realized when σ ≤ κ and a growth solution when σ > κ .The value of Vi at the growth solution is obtained using the nonlinear equation, which is obtained by eliminating θ from Eqs. (3) and (4).
The solution to this equation, Vi , is obtained by numerically searching between V * and V∞ using the bisection scheme.

Diffusion of solute and impurity molecules in solution
We assumed that the crystal-solution interface is an infinitely wide plane and is present on the x-axis from the interface perpendicularly toward the solution.The molar concentration distribution of the solute and impurity molecules in the solution at a certain time t is denoted by c α (x, t) , where α denotes the type of molecule, and the solute and impurity are denoted by α = s and α = i , respectively.The solute and impurity are depleted near the interface by being incorporated into the growing crystal and then supplied by diffusion from the far side.The where D α is the diffusion coefficient.We assumed that the crystal growth is sufficiently slow and that the move- ment of the interface does not affect c α .In addition, when solving Eq. ( 11), the position of the interface is always assumed to be x = 0.
The boundary condition at the interface is determined by the conservation law for each molecule.The number of moles of molecules flowing into the interface from the solution per unit area and unit time (interfacial flux) is given by Fick's law of diffusion: j α = D α ∂c α ∂x .If the impurity fraction with respect to the solute is sufficiently small, we can assume that the crystal thickens by the amount of solute incorporated into the interface.The relationship between the normal growth rate R (the crystal growth rate perpendicular to the interface) and interfacial flux of solute j s is expressed as R = v m j s .Thus, we can obtain the following boundary condition for solute molecules: Conversely, the amount of impurity molecules incorporated into the crystal is determined using the partition coefficient K, which is the ratio of the fraction in the crystal to that in the solution.The impurity fraction in the solution near the interface is expressed as c i (0, t)/c s (0, t) , and that in the grown crystal is expressed as j i /j s .Therefore, the distribution coefficient is expressed as K = j i /j s c i (0,t)/c s (0,t) .Thus, we can obtain the following boundary condition for impurity molecules: In this study, we considered a closed system.The size of the computational domain is set to L and the no-flux boundary condition ∂c α ∂x = 0 is imposed at x = L .L was set to be large enough to ensure that the effect of the boundary layer does not reach it.

Calculation procedure
The interfacial supersaturation was calculated as σ = c s (0,t)−c e c e , where c e is the solubility.As two solutions of V exist when σ * < σ < κ , we used a logical variable w, which represents the growth state of the crystal surface, to determine which solution is chosen.The initial value of w was set as "True;" changed to "False" when σ fell below σ * and to "True" when σ exceeded κ .When w was "True, " we solve Eq. ( 10) to obtain Vi , and when it was "False, " we just set V = 0 .The normal growth rate was determined using R = (a/ )V i = (a/τ ad )σ Vi , where a is the step height.Note that the calculation result does not depend on the value of 0 .R was substituted into the boundary conditions (Eqs.12 and 13) and the diffusion equation (11) was solved numerically to obtain the time variation of the concentration distributions c α (x, t) .The computational domain from x = 0 to L was divided by cells with equal intervals of x , and Eq. ( 11) was solved using the finite volume method.The time increment was set to t , and the first-order upwind difference was used for the time derivative.By repeating these procedures, we obtained the time variations of R(t) and c α (x, t) .The growth thickness is given by h(t) = t 0 R(t ′ )dt ′ and the impurity fraction of the grown crystal surface is given by χ(t) ≡ j i j s = K c i (0,t) c s (0,t) , so the impurity distribution in the grown crystal is shown in the plot of h(t) versus χ(t).
The values of the parameters used in the calculations are summarized in Supplementary Table 1.Herein, typical solution growth situations were considered.The diffusion coefficients of the solute molecules were taken from typical values for metallic ions in aqueous solution 31 .Impurity molecules were assumed to be larger than solute molecules and their diffusion was assumed to be slower than that of solute molecules.For solubility, typical values for soluble substances in water were used 32 .The initial solute concentration was set sufficiently high to allow crystal growth.The initial concentration of impurity molecules was 1/100 of the solute molecules.The adsorption time τ ad of the impurity molecules on the growing crystal surface is not well-known.The adsorp- tion time of macromolecules, such as proteins, on inorganic crystal surfaces varies from less than one minute to several tens of minutes, depending on the inorganic crystal and impurity molecules 16 .Analyzing ice crystal growth hysteresis with synthetic antifreeze proteins has revealed extremely short adsorption times of < 0.1 s 21 .With reference to the literature, the case of τ ad = 0.1 s is considered in this study.For the inorganic crystals, the step kinetic coefficient is estimated as β st ∼ 10 −2 -10 −1 cm/s 30 ; thus, V 0 = β st v m c e ∼ 10 −5 m/s .Assuming that the step interval is 0 ∼ 100 nm , V0 = V 0 /( 0 /τ ad ) ∼ 10 .Thus, we considered the case where V0 = 5 .Assuming that the halted growth of a crystal will resume when σ > 1 , we set κ = 1 .L was sufficiently large to discontinue the influence of the boundary layer on the outer edge of the computational domain.Our findings for the slow impurity adsorption ( τ ad = 0.2 s and V0 = 10 ) case are shown in Fig. 3, confirming its dependence on impurity adsorption kinetics.

Width of the bistability region
The supersaturation σ * at the transition point is dependent on the value of V0 κ2 (Eqs.8 and 9).For V0 κ2 < 1 , σ * / κ is almost 1 20 , indicating that the bistability region becomes very narrow, and therefore growth hysteresis is unlikely to occur.Here we have The cases where the value of V0 κ2 is small are as follows.First, c e is small (insoluble), τ ad is small, or β st is small.This corresponds to a situation where crystal growth is slower than impurity adsorption, and adsorption/desorption equilibrium is always reached at the crystal surface, thus, the interdependence between step velocity and adsorbed impurity density is untenable.Second, N e is small and step advancement cannot be halted at inter- mediate supersaturation conditions because the adsorbed impurity density is low at the adsorption/desorption equilibrium.Hence, the bistability region only appears under limited supersaturation conditions; thus, growth hysteresis is unlikely to occur, suggesting that no OCZ is generated by this mechanism.

Figure 1 .
Figure 1.Mechanism of growth inhibition by impurities adsorbed on the crystal surface.(A) Interaction between adsorbed impurities and atomic steps.The crystal surface has a step, which is as high as one atom, and the continuous incorporation of atoms into this step moves the step forward and builds the surface layer by layer (layer growth).The step is generated by screw dislocations exposed to the surface.Adsorbed impurities inhibit the step advancement by pinning it in place (pinning effect).(B) Steady-state solution of the step velocity when considering impurities that are repeatedly adsorbed and desorbed on the surface.The abscissa and ordinate represent the degree of supersaturation and dimensionless step velocity, respectively.In an intermediate supersaturation region (bistability region), there exist two solutions: growth (purple) and no-growth (blue) states.The dashed curve indicates an unstable solution that does not appear in practice.

Figure 2 .
Figure 2. Calculated normal growth rate and concentration distributions in solution.(A) Time variation of the normal growth rate.(B) Enlargement of abscissa of panel A. (C) Time variation of concentration distributions in the solution during crystal growth.The solid and dashed lines indicate the solute and impurity concentration distributions, respectively.The impurity concentration is multiplied by a factor of 100 for comparison with the solute concentrations.Labels a-c correspond to those in panel B. (D) Time variation of the concentration distributions when the crystal stops growing.Labels d-f correspond to those in panel B. In this calculation, impurity diffusion is assumed to be slower than the solute diffusion.

Figure 3 .
Figure 3. Impurity distribution in the grown crystal.The abscissa represents the distance from the initial crystal surface position, and the ordinate represents the impurity fraction relative to solute molecules.The blue line shows the case in Fig.2.The orange line shows the case when the impurity adsorption is assumed to be slower.
σ )convective transport is ignored.The time variation of c α is expressed by the following unsteady one-dimensional diffusion equation: